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Motivated by recent efforts to achieve cold fermions pairing near a Feshbacli resonance, we 
consider tfie dynamics of formation of tfie Bardeen-Cooper-Schrieffer (BCS) state. At times shorter 
than the quasiparticle energy relaxation time, after the interaction is turned on, the dynamics of the 
system is nondissipative. We show that this collective nonlinear evolution of the BCS-Bogoliubov 
amplitudes Up, Vp, along with the pairing function A, is an integrable dynamical problem, and 
obtain a family of exact solutions in the form of single solitons and soliton trains. We interpret 
the collective oscillations as Bloch precession of Anderson pseudospins, where each soliton causes 
a pseudospin 2tt Rabi rotation. Numerical simulations demonstrate robustness of the solitons with 
respect to noise and damping. 



Dilute fermionic alkali gases cooled below degener- 
acy [1] are expected to host the paired BCS state [2,3]. 
One of the unique and attractive features of this system 
which makes it qualitatively different from superconduct- 
ing metals is the possibility to control the strength of 
pairing interaction and change it by using the Feshbach 
resonances [4-9]. Pairing enhancement near the reso- 
nances, as well as the high coherence of atomic systems, 
can allow to realize the strong coupling BCS regime, as 
well as to explore the time dynamics of the paired state. 
Since the characteristic energy scales in atomic vapors 
are relatively low, one can perform time-resolved mea- 
surements on the intrinsic microscopic time scales, and 
explore a range of fundamentally new phenomena in the 
time dynamics of the paired state. These new prospects 
helped to revive interest in some of the basic issues of the 
BCS pairing problem [10-15]. 

In particular, one of the important questions has to 
do with the dynamics and characteristic times of forma- 
tion of the BCS state in cold gases [16]. The dynam- 
ics of the superconducting state in metals, described by 
BCS theory, has been a subject of active research for a 
long time [17]. Generally, there are two important time 
scales to be considered: the quasiparticle energy relax- 
ation time Te, and the characteristic time of the order pa- 
rameter change, ta, estimated as the inverse increment 
of Cooper instability [18,19]. For <C ta, quasiparti- 
cles quickly reach local equilibrium parameterized by a 
time-dependent order parameter A{t). In this case the 
dynamics is described by the time-dependent Ginzburg- 
Landau (TDGL) equation for A{t), with the relaxation 
time scale ca. ta. However, as noted by Gorkov and 
Eliashberg [20], this scenario applies only to relatively 
exotic situations, including a close proximity of a transi- 
tion point, A^/Tc <^ T~^, or a fast pair breaking (e.g., 
due to paramagnetic impurities). 

In the opposite limit, which takes place at the temper- 
atures not too close to critical, 

Te > TA (1) 

the system dynamics is usually described by Boltzmann 
kinetic equation for quasiparticles and a self-consistent 



equation for A(t) connecting it with the quasiparticle 
distribution function locally in time [21,22]. The valid- 
ity of this approach requires, in addition to (1), that the 
variation in time of both the quasiparticle distribution 
and the external parameters is sufficiently slow on the ta 
time scale. The origin of this criterion is seen most easily 
if one notes that the order parameter can exhibit free os- 
cillations with a frequency ca. [23] (also, see below). 
Thus, at a slow parameter variation the order parame- 
ter A(t) and the quasiparticle spectrum will adiabatically 
follow the changes of the quasiparticle distribution, with- 
out exciting the oscillations of A. This approximation is 
relevant to the majority of situations in superconducting 
metals because of the relatively big value of A and the 
difficulty of changing the external parameters and mak- 
ing measurements on the time scale of ta. 

From this point of view, the cold fermionic gases 
present a completely new situation. The energy relax- 
ation in these systems is quite slow, while the external 
parameters, such as the detuning from Feshbach reso- 
nance, can change very quickly on the time scale of ta. 
Thus the BCS correlations in this case build up in a co- 
herent fashion while the system is out of thermal equilib- 
rium. In such a situation, theory must account not only 
for the order parameter evolution, but also for the full 
dynamics of individual Cooper pairs and quasiparticles. 
It seems to be desirable to understand better this 'fast 
BCS buildup' regime, since the lifetime of the gas sam- 
ples is finite, while the time-resolved measurements can 
easily be performed with resolution better than ta. 

In this article we consider the situation when the pair- 
ing interaction is switched on essentially instantly, on a 
time scale tq <C ta , t^ . We shall discuss only the spatially 
uniform situation, relevant for samples of finite size, and 
explore the time evolution of the unpaired ideal Fermi 
gas, which is unstable with respect to Cooper pairing. At 
not too long times, t ^ t^, the dynamics of the system 
is governed by nondissipative equations which conserve 
both the entropy and the total energy. A stationary su- 
perconducting state with the same energy as that of the 
initial state would have more quasiparticles, and thus 
would have entropy higher than the initial Fermi gas en- 
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tropy. This argument suggests that the system can reach 
a stationary (not necessarily equiUbrium) state only at a 
very long time of the order t^. On shorter times, which 
nevertheless can exceed ta), the system will exhibit a 
nonlinear time evolution. During this period of time the 
concept of the quasiparticlc spectrum is irrelevant and 
theory can rely neither on the kinetic equation, nor on 
TDGL equation. Below, we present an approach which 
describes the BCS state buildup and accoimts for coher- 
ent dynamics of individual Cooper pairs. We shall focus 
on the zero temperature case, when ta — A~^, and show 
that the result of Cooper instability is a periodic oscilla- 
tion A(t) having the form of a soliton train. 

This regime can be described by the BCS hamiltonian 



X{t) 



pi^-qi'^qT ' 



(2) 



with the coupling turned on abruptly, X{t) = X9(t — i,). 

The main result of this work is that the time-dependent 
problem (2) is integrable. We generalize the BCS solu- 
tion, which is exact for the separable pairing Hamiltonian 
(2), and demonstrate that at f > the many-body state 
evolves as 



= n («p(*) +%(*)<T«^p,i) io) 



(3) 



The Bogoliubov mean field treatment, which gives a state 
of the form (3), relies on the 'infinite range' form of the 
pairing interaction in (2) (i.e., equal coupling strength of 
all (p, — p), (q, — q)). Since the latter does not depend 
on the system being in equilibrium, one can introduce a 
time-dependent mean field pairing function 



A{t) = xJ2^^itKit) 



(4) 



The amplitudes Up(t), Vp{t) can be obtained from the 
Bogoliubov-deGennes equation 



idt 



ep A 

A* -e. 



(5) 



to be solved along with the selfconsistency condition (4). 

We recall that the unpaired state is a selfconsistent, al- 
beit unstable, solution of Eqs. (5), (4) with A = 0, T = 0: 

u^^\t) = e-^^p*6'(ep), v^\t) = e*^p*6l(-ep). The stability 
analysis [19] shows that the deviation from the unpaired 
state grows as A{t) oc e'''*e~*'^*, 



6up{t) = 



A{t)v^°\t) 



(6) 



ij - 2ep -h w Z7 -h 2ep - oj 

with the growth exponent 7 and the constant u) given by 
sgncp 



2ep-C' 



( = U! + ij 



(7) 



The electron-hole symmetry near the Fermi level, 
N[ep) = N{—ep), makes the frequency shift w vanish. 
Using the similarity between Eq. (7) and the BCS gap 
equation at T = 0, the exponent 7 can be shown to be 
close to the BCS gap value Aq. (In fact, in the constant 
density of states approximation, the two quantities coin- 
cide: 7 = Aq.) Thus, we estimate the time ta — Aq ^. 

The interaction switching, while nonadiabatic, must 
also be gentle enough not to overheat the fermions above 
T ~ Aq. The effective temperature T^g after switching 
can be estimated from the total energy increase, giving 

oT^ff = Y^hojnin2{l - n3)(l - 714)1X^1"^ S{hoj - SE) (8) 



w, 1...4 



,, a 



^v. For 
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with rii = n(epj, 5E = Cpg+ep^-epi 
the interaction switching from to A over a characteristic 
time To, the RHS of (8) is of the order of A^Tq"^. Thus the 
cold Teff condition is EfTq ^ (An/Ag)^/^, which is com- 
patible with the nonadiabaticity requirement tq ^ ta- 

At T = 0, a soliton solution of Eqs. (5), (4) can be con- 
structed most naturally in terms of the variable 



w ={ "p/^P' ^P > ^ 
P 1 Vp/up, Cp <0 



(9) 



Consider first the case Cp > 0. Prom Eq. (5) we obtain 

idtwp = 2epWp + A{t) - A* {t)wl (10) 

where A{t) is a function of time determined from (4). 

Motivated by the stability analysis (6), (7), we try the 
ansatz A{t) = e~^'^^a~^{t) with a{t) real, and 

Wpit)=2epf{t)-zf{t), /(i) ^ i- = e-^'^*a(t) (11) 

Substituting this in Eq.(lO), we obtain an equation 

i^pct + a = ^p(Cpa -ia) + ^{l- {^pU - iaf) (12) 

with ^p = 2ep— a;. Remarkably, the terms with ^p cancel, 
and Eq. (12) takes the same form for all the states, 



aa = d +1 



(13) 



which justifies the ansatz (11). By a variable substitution 
a = e'^, Eq. (13) can be brought to the form (p = e~^'^. 

Integrating the latter equation, obtain (f)'^ + e~^* = 7^, 
with 7 an integration constant. This yields d^ = 7^a^ — 1, 

a(t) = -cosh7(^-to), m = ^J^^,_,^^ (14) 

The modulus |A| time dependence (Fig.l, upper left), 
growing first, then decreasing and taking the system back 
to the unpaired state, reflects the absence of dissipation. 
The time to at which |A| peaks is set by the initial con- 
dition at large negative times t ~ . 
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For ep < 0, Wp = Vp/up, the form of Eq. (10) remains 
the same up to a sign change Cp — > —Cp and the permuta- 
tion A ^ A*. Accordingly, the ansatz for Wp in this case 

is w.^Mt) = 2|ep|/W - f{t) = A-i = e'-*a(i). 

The resulting equation for a{t) is identical to Eq. (13). 

The last step is to analyze the requirements on this so- 
lution due to the selfconsistency condition (4). For that, 
we rewrite Eq. (4) in terms of Wp{t) as 



Wp{t) 



ep>0 



+ A 



Cp<0 



(15) 



and note that both the right and the left hand side have 
the same time dependence, and are equal to each other 
provided that the quantity C, = u + satisfies Eq. (7). 
This means that the gap modulus |A(t)| peak value is 
equal to the instability exponent 7 defined by (7). 
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FIG. 1. Coherent BCS dynamics on the Bloch sphere (16). 
Trajectories of individual Cooper pair states for the soliton 
train solutions (Eq. (26) and insets) oscillating in the limits 
A- < A{t) < A+ are shown. Note that each state completes 
a full 2n Rabi cycle per soliton. The red and blue curves cor- 
respond to the energies ep above and below the Fermi level. 

At r = 0, in the case of a constant density of states, 
7 is equal to the equilibrium BCS gap Aq, while uj van- 
ishes due to particle-hole symmetry. Thus, remarkably, 
the modulus |A(t)| in this case peaks exactly at Aq. 

To illustrate the collective dynamics in the soliton so- 
lution, we plot the trajectories Up{t), Vp{t) on the Bloch 



sphere 



+ — 1 using the parameterization 



ri + ir2 



2upVp; 



(16) 



(Fig. 1, left). Here, as well as in the soliton train solutions 
discussed below (Eq. (26)), each state {up,Vp) completes 
a full Rabi cycle per soliton. The trajectories, which are 
small loops for the pairs with large energies Cp, turn into 
a big circle as ep tends to the Fermi level. 

To gain further insight, we reformulate the Bogoli- 
ubov approach, following Anderson [23], in terms of pseu- 
dospins associated with individual Cooper pair states. 
'Pauli spin' operators = ^(Cp iicr^) can be assigned 
to each pair of fermions with opposite momenta as follows 



flp^a. 



Pi 



-p1«pT 



(17) 



allows 



The identity = [cr+,crp] = a'^^ap-^ - a^p^c^^^ 
to map the problem (2) onto an interacting spin problem 



p p,q 



(18) 



where means a sum over the pairs of states (p, — p). 
Since all the spins interact with each other equally, the 
mean field theory here is exact, just like for the BCS 
problem. The mean field hamiltonian for each spin is 



(-A',-A",6, 



(19) 



Here the z component of the effective field bp, given by 
the single particle energy, is spin-specific, while the trans- 
verse components, the same for all of the spins, satisfy 



(20) 



which is analogous to the BCS gap relation. In the 
ground state each spin is aligned with bp, and the spins 
form a texture near the Fermi surface [23] , with spin ro- 
tation described by the Bogohubov angle. 

The dynamical problem of interest can be cast in the 
form of Bloch equations for the spins. 



(Tp = «[7Yp,(Tp] = 2bp 



(21) 



with the field bp defined selfconsistently by (19), (20). 
Anderson [23] used Eq. (21), linearized about the texture 
ground state, to study collective excitations of a super- 
conductor (see also [24]). Linearized about the unpaired 
state, Eq. (21) yields an instability identical to (6), (7). 

The Bloch dynamics is suitable for simulation (Fig. 2), 
since Eq. (21) is linear in spin operators and can be re- 
placed by an equation for the expectation values (16). 
Typical A{t), observed in the presence of thermal noise, 
is an orderly sequence of the cosh solitons (Fig. 2, top), 
indicating robustness of soliton solution (c/. Ref. [25]). 
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FIG. 2. Bloch dynamics Vp — 2bp x rp for 10 spins with a 
constant density of states, with bp*^ defined by (19), (22). Ini- 
tial conditions r^^p = tanh(i/3ep), {ri+ir2)p = e"''p(l — rf) 
with /3~^ = T — O.lAo and random uncorrelated < 
were used to simulate thermal noise. 
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The effect of damping can be studied by replacing 



bp = bp - 7bp X Tp 



(22) 



with a dimensionless damping constant, 7 ~ l/(AoTe). 

This destroys intcgrability and prompts relaxation to 
BCS equilibrium (Fig. 2, bottom). 

In our nonlinear dynamical problem, we have to solve 
the Bloch equations (21) together with the selfconsis- 
tency relation (20). As above, wc assume the pairing am- 
plitude time dependence of the form A{t) = e~^'^*il{t), 
where Q{t) is real and w is a constant frequency shift to 
be fixed by the relation (20). In the Larmor frame rotat- 
ing with the frequency lu, Eq. (21), written for the spins 
average polarization = {a^} components, reads 



n 



<pr2, r2 = ^pri + 2fir3, rs = -20r2 (23) 



with ^p = 2ep — u), as before. An exact solution can be 
obtained from the ansatz 



n = An, r2 = Bfl, rz = CCl^ - D 



(24) 



The first and the third equation (23) are satisfied by (24) 
provided A = —^pB and B = — C, while the second equa- 
tion (23) is consistent with the normalization condition 
rf -|- -|- = 1, and thus yields 

C^C^O^ + C^tf + (CQ2 -Df = l (25) 

Eq. (25) will take the same form for all the spins, 

J72^(02-A?_)(J72_A^)=0, A_<A+ (26) 

provided that the constants D, C are chosen as 

(£>2-l)/C2 = A^ A^, 2£>/C = Cp + A?. A^ (27) 

with the sign factor sgnep. Eq. (26) defines an elliptic 
function rj(t) oscillating periodically between A_ and 
A+. At A_ <C A+, the solution is a train of weakly 
overlapping solitons (14) with A+ = 7 (Fig. 1, left). 
The real part of the selfconsistency relation (20), 



CpSgnep 



((^2+A2 +A^)2-4A2A2) 



1/2 



(28) 



fixes one of the constants A±, the other one being a free 
parameter that controls the inter-soliton time separation. 
With the ratio A_/A_|_ varying from to 1, the soliton 
frequency increases, and the weakly overlapping solitons 
(14) gradually become overlapping more strongly, turn- 
ing into weak harmonic oscillations (Fig. 1). 

We note that the weak oscillation limit has been inves- 
tigated by Volkov and Kogan [24] with the help of lin- 
earized Gorkov equation, and earlier by Anderson [23], 
who used pscudospin representation. The nonexponen- 
tial decay of the lianearized oscillations [18,24] was inter- 
preted as coUisionless damping, related to strong mixing 



of the oscillations of A with the states of two excited 
quasiparticles slightly above the superconducting gap. 

The imaginary part of Eq. (20) fixes the value of the 
frequency shift w (we recall that w 7^ in the presence of 
charge asymmetry). At A_ <C A+, Eq. (28) turns into 
Eq. (7) which, as we found above, defines the amplitude 
of a single soliton. In the opposite limit, A_ — »■ A+, 
Eq. (28) coincides with the BCS gap equation. 

There is an interesting relation between our problem 
and the self- induced transparency phenomenon [26]. In 
the latter, an optical pulse interacting with an ensemble 
of atoms can dissipate its energy by inducing resonant 
Rabi transitions in the atoms. However, when the pulse 
duration is tuned so that the atoms complete a full Rabi 
2ti cycle as the pulse goes by, the pulse sustains itself and 
propagates without dissipation. The Bloch equations de- 
scribing this phenomenon bear striking similarity to our 
Eqs. (21), while the atoms' polarization is employed [26] 
to provide feedback on the pulse instead of our BCS self- 
consistency relation. 

Before concluding, we note that the dynamics at finite 
temperature, in the regime described by Eq. (1), remains 
an open problem. In particular, we can not rule out 
the possibility of chaotic behavior. At T = the prob- 
lem has a relatively simple solution, periodic in time, be- 
cause in this case the quasiparticles with low energies are 
strongly coupled to the oscillations of A. In contrast, in 
the case Tta ^ 1 there are two groups of quasiparticles: 



those with energies ep 



which fully participate 



in the oscillations, as above, and the quasiparticles with 
<C Ep <C T, coupled to A(f) much weaker and playing 
the role of a thermal bath, thereby providing dissipation. 

In summary, this work provides an exact solution for 
the BCS pairing formation problem. In the nonadiabatic 
regime, the dynamics is dissipationless and thus nonlin- 
ear. Soliton train solutions are obtained analytically and 
demonstrated to be generic and robust by a simulation. 
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